clc
clear

stencil_width = 3;

syms dx
q = sym('q',[1 stencil_width]);

qL = WENO3(q,-1,dx);
qR = WENO3(q, 1,dx);

for iCell = 1:stencil_width
    dqLdq(iCell) = diff(qL,q(iCell));
    dqRdq(iCell) = diff(qR,q(iCell));
end

ccode([dqLdq;dqRdq],'file','WENO3_jab.c');

function qrec = WENO3(q,dir,dx)
syms epsilon
C = sym('C',[1 2]);

if(dir>0)
    qim(1) = q(1);
    qi     = q(2);
    qip(1) = q(3);
else
    qim(1) = q(3);
    qi     = q(2);
    qip(1) = q(1);
end

stencil(1) = -1./2. * qim(1) + 3./2. * qi;
stencil(2) =                   1./2. * qi + 1./2. * qip(1);

coefA(1) =-qim(1) + qi;
coefA(2) = qip(1) - qi;

beta = coefA.^2;

% WENO-Z
tau  = abs( beta(2) - beta(1) );
xi = tau ./ ( beta + epsilon );
alpha = C .* ( 1 + xi.^2 ); % WENO-Z
%alpha = C .* ( 1 + tau ./ ( beta + epsilon ) + dx^(2/3) * beta / ( tau + epsilon ) ); % WENO-Z+

% WENO-JS
% alpha = C ./ ( beta + epsilon ).^2;

omega = alpha ./ sum(alpha);

qrec = sum(stencil.*omega);

end
